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A non Monte Carlo path-integral algorithm that is particularly adept at handling nonlinear 
Lagrangians is extended to multivariate systems. This algorithm is particularly accurate for systems with 
moderate noise. 
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PATHINT is a non Monte Carlo histogram C computer code developed to evolve an n -dimensional 
system (subject to machine constraints) based on a generalization of an algorithm demonstrated by 
Wehner and Wolfer to be extremely robust for nonlinear Lagrangians with moderate noise [1-3]. 
PATHINT was recently used in a neuroscience study [4], where it was observed how difficult it can be 
even in two dimensions to get good resolution because of CPU constraints on Sun SPARC 2 and 10MP 
machines. Here it appears that the resolution is quite satisfactory on these machines. 

The system selected for this paper to illustrate the use of PATHINT is the classical analog of a 
quantum system. "Quantum chaos" was a term used to describe the observation of chaos in the classical 
trajectories of the Hamilton's equations of motion, 

• -_ d JL 
^ dq 



dH 

4=^, (1) 
op 

where a Hamiltonian of the form 
H= 1 -p 2 +4>(q), 



<J> = - (VO) 2 - - DV 2 <I> , (2) 

2 2 

with a coefficient D was considered in a q dimension > 2, representing double this dimension in the 
corresponding phase space of q and p [5]. 

In ref. [5] a potential was considered, 

O = 2x 4 + ^ y 4 + exy(x - yf , (3) 
in a classical Fokker-Planck system 

dP, = V. (PV<D) + y V 2 P . (4) 
A slight change of notation is used in this paper. They posed the question as to just what properties such 
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classical systems might possess? 

This paper does not at all deal with the quantum system described above, but it does deal in detail 
with the associated classical system, and it does give some answers to the above posed question. This 
study should be considered as illustrating a particular numerical approach that promises to be quite useful 
in studying the evolution of such classical systems. The results obtained can be considered as 
"experimental" data on the exact region of such classical transformations. The results here are negative 
with respect to any unusual or interesting activity in the parameter region observed in the quantum 
mechanical calculation. This should at least help other investigators who might tend to focus on this 
region in the classical system, based on the results obtained for the associated quantum system. 

This paper computes the path integral of the classical system in terms of its Lagrangian L. 

( ' ) 

P[q t \q to ]dq(t) = j ■■■ j Dqexp -mm jdt'L S(q(t ) = q )S(q(t) = q t ) , 

V *o J 



H+1 

Dq= lim ng 1/2 ri(2^A0- 1/2 ^ ! p , 



L(q i ,q i ,t)= 1 -(q i -g i )g ir (q' ~ /) , 



g = det(g, r ) . (5) 

Here the diagonal diffusion terms are g xx = g yy = D and the drift terms are g l =-d<frldq'. If the 
diffusions terms are not constant, then there are additional terms [6]. 

The histogram procedure recognizes that the distribution can be numerically approximated to a high 
degree of accuracy by sums of rectangles of height P ; and width Aq' at points q'. For convenience, just 
consider a one-dimensional system. The above path-integral representation can be rewritten, for each of 
its intermediate integrals, as 
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P(x; t + At) = \ dx'[g U2 (2xAty m exp(-L At)]P(x'; t) 
= J dx'G(x, x'\ At)P{x'\ t) , 



P(x; t) = £ x{x - x ! ')P,<0 , 
1=1 



n(x-x l ) 



1 . (x ! - - Ax' -1 ) < x < (jc 1 + - Ax 1 ) , 

2 2 (6) 

0, otherwise. 



This yields 

P i (t + At) = T ij (At)Pj(t), 

TJAt) = — r-j r f dx f ^'G(x, jc'; A?) . (7) 

J Ax'" 1 + Ax ! J x'-Ax^n J xJ-AxJ- l /2 

Tjj is a banded matrix representing the Gaussian nature of the short-time probability centered about the 
(possibly time-dependent) drift. 

This histogram procedure was extended to two dimensions using a matrix 7^ [7]. Explicit 
dependence of L on time t also can be included without complications. Care must be used in developing 
the mesh in Aq l , which is strongly dependent on the diagonal elements of the diffusion matrix, e.g., 

Aq l - {Atg U ) m . (8) 

This constrains the dependence of the covariance of each variable to be a (nonlinear) function of that 
variable in order to present a straightforward rectangular underlying mesh. 

Since integration is inherently a smoothing process [8], fitting the data with the short-time 
probability distribution, effectively using an integral over this epoch, permits the use of coarser meshes 
than the corresponding stochastic differential equation. For example, the coarser resolution is 
appropriate, as typically required, for a numerical solution of the time-dependent path integral. By 
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considering the contributions to the first and second moments, conditions on the time and variable meshes 
can be derived [1]. The time slice essentially is determined by At < L l , where L is the uniform 
Lagrangian, respecting ranges giving the most important contributions to the probability distribution P. 
Thus At is roughly measured by the diffusion divided by the square of the drift. 

Such calculations are useful in many disciplines, e.g., some financial instruments [8,9]. Monte 
Carlo algorithms for path integrals are well known to have extreme difficulty in evolving nonlinear 
systems with multiple optima [10], but this algorithm does very well on such systems. The PATHINT 
code was tested against the test problems given in previous one-dimensional systems [1,2], and it was 
established that the method of images for both Dirichlet and Neumann boundary conditions is as accurate 
as the boundary element methods for the systems investigated. Two-dimensional runs were tested by 
using cross products of one-dimensional examples whose analytic solutions are known. 

Attempts were made to process the same system considered for the quantum case [5]. Therefore, 
the diffusion was taken to be D = 0. 2. Since they selected an harmonic oscillator basis for their 
eigenvalue study, is was assumed that natural boundary conditions are appropriate for this study, and 
ranges of x and y were tested to ensure that this was reasonable. A band of three units on each side of the 
short-time distribution was sufficient for these runs. For e < 5, the range of x was taken to be +3 and the 
range of y was taken to be ±6. 

A mesh of At = 0. 1 was reasonable to calculate the evolution of this system for 0. 1 < e < 0. 5. The 
quantum study observed chaos at e > 0. 14, but the classical system appears to be very stable with a single 
peak in its probability density up through s = 0. 5. The time mesh was tested by performing several 
calculations at time meshes of At = 0. 01 on a Sun SPARC 10MP. All other calculations reported here 
were performed on a Sun SPARC 2. 

Figures 1(a) and 1(b) show the evolution of the distribution for e = 0.1, t = 10 and 100, i.e., after 
100 and 1000 foldings of the path integral, respectively. The distribution starts at single peaks at the 
imposed initial condition, a 8 function at the origin, i.e., x = -0. 0302 and y = -0. 0603 with this mesh, 
and swell out their stable structures within t < 0. 03. Note the stability over the duration of the 
calculation. A similar stability was noted for e < 0.55, as illustrated in Figures 2(a) and 2(b) for e = 0.3 
and 0.5 after 1000 foldings. 
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For values of e > 0. 5, a mesh of At = 0. 05 was used to investigate the onset of instabilities which 
were noted for higher e with coarser meshes. To two significant figures, this first occurs at e = 0. 55. As 
e increases, so does the peak spread quite early along a diagonal in the evolving distribution. As the time 
of the calculation increased, there was concern that the boundaries were being approached. Therefore, for 
these runs, the range of x was increased to ±5 and the range of y was increased to +10. Figure 3 
illustrates the evolution at time t = 5. (100 foldings) and t = 100 (2000 foldings) with e = 0.55. Figure 4 
illustrates the early structure for e = 0.6 at t = 15 after 300 foldings of At = 0. 05. 

The PATHINT algorithm utilized here will be used to explore the time evolution of other Fokker- 
Planck systems in the presence of moderate noise. As mentioned above, such problems arise in many 
systems ranging from neocortical interactions to financial markets. Also, we can now accurately examine 
long-time correlations of chaotic models as multiplicative noise is increased to moderate and strong 
levels; many chaotic models do not include such levels of noise as is found in the systems they are 
attempting to model. A project is now underway under an award of Cray computer time from the 
Pittsburgh Supercomputing Center through the National Science Foundation (NSF), the Parallelizing ASA 
and PATHINT Project (PAPP) further developing this code for large systems [11]. 
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FIGURE CAPTIONS 

FIG. 1. Probability density for e = 0. 1, for At = 0.1 (a) after 100 foldings (t = 10) and (b) after 
1000 foldings (t = 100). 

FIG. 2. Probability density after 1000 foldings of At = 0A(t = 100) for e set to (a) 0.3 and (b) 0.5. 

FIG. 3. Probability density for e set to 0.55 using At = 0. 05, (a) at t = 5 after 100 foldings, and (b) 
at t = 100 after 2000 foldings. As discussed in the text, the ranges of x and y for s > 0. 5 were increased. 

FIG. 4. Probability density at t = 15, after 300 foldings of At = 0. 05, for e = 0.6. 
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Figure 2a 
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Figure 3a 



Pith- integral evo I utio n of m univariate syEttniE Lester Lnjber 

Flpure3j3 



i I 




cor-uaLn^cicy^-o 
dddddddd 



Path-integral evolution of multivariate systems 



Lester In; 



Figure 3b 

Pith- integral evo I utio n of m univariate syEttniE Lester Lnjber 

Figure 3b 




Path-integral evolution of multivariate systems 



Lester In; 



Figure 4 
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